function [results, SBEdata] = FlightsDriver_dev(varargin)
  % Driver for MMU flights data processing
  %
  % Author: David Goldsmith, Wash. State Dept. of Ecology, dgol461@ecy.wa.gov
  % Original release date: 9/26/2011

    global lfcr P T C dOV PH xmiss Chl NTU par S dOS DO v 
    
    if isempty(lfcr)
        lfcr = char(10);
    end
    
    if ~is_devVer(mfilename)
      % Usage
        input(['Welcome to FlightsDriver (FD), a driver for MISU''s re-engineering of Seabird''s' lfcr...
               'SBE Data Processing (SBEDP) software package! Presently, FD reads an SBE25 raw' lfcr...
               'hex file and an SBEDP xmlcon file of the user''s choosing, and produces MATLAB' lfcr...
               '"vectors" (i.e., 1-D arrays) of the recorded pressure, temperature, conductiv-' lfcr...
               'ity, pH, transmission, fluorescence/chlorophyl, turbidity, and PAR, as well as' lfcr...
               'the calculated salinity, dissolved oxygen saturation concentration, and dissol-' lfcr...
               'ved oxygen voltage. Then it offers the user the option of graphically comparing' lfcr...
               'its results to those of SBEDP (produced independently but using the same input' lfcr...
               'hex and con files, of course). (Note: if run from within MATLAB, both FD''s and,' lfcr...
               'if you chose to perform the above comparison, SBEDP''s results can be returned to' lfcr...
               'the command window or a calling function, but the compiled binary version exits' lfcr...
               'after producing the comparison plots.)' lfcr lfcr...
               'Now kindly press Enter to continue...' lfcr]);
    end
    results = {}; files = {}; paths = {};
    if nargin==0 % Use a file chooser dialog to get required hex and xmlcon files
        while numel(files) < 2
            [filesTmp, pathsTmp, ~] = uigetfile('*.hex; *.xmlcon',... % Launch file chooser dialog
                                                 'Hex/xmlcon file chooser',...
                                                 'MultiSelect', 'on');
            if ~iscell(filesTmp), filesTmp = {filesTmp}; end % If only a single file is chosen, 
                                                             % put it in a cellarray for API
                                                             % uniformity
            files = extendArray(files, filesTmp, 'c');
            for i=1:numel(filesTmp)
                paths{end+1} = pathsTmp;
            end
        end
    else
      % TODO: Implement command-line input of file names  
    end
    
  % Create full paths to the required files
    idx = ~cellfun(@isempty, regexpi(files, '.hex$'));
    hexFile = fullfile(paths{idx}, files{idx});
    idx = ~cellfun(@isempty, regexpi(files, '.xmlcon$'));
    conFile = fullfile(paths{idx}, files{idx});

  % Instantiate data container
    results = ParseRawSBE25Data_dev(hexFile, conFile);
    
  % Parse headers for instrument data
    disp('Parsing headers')
    results.parseInstHeader;
    
  % "Unpack" configuration data
    disp('Unpacking configuration data')
    results.unpackConfigDat;
    
  % Convert hex sensor data to decimal
    disp('Converting hex to decimal')
    results.convertHexData;

  % Convert decimal sensor data to physical units
    disp('Converting decimal raw sensor data to "physical" units')
    disp('Opening connection to EAPMW')
    conn = database('EAPMW', '', '');
    disp('Computing pressure')
    P = Pressure_dev(results, conn); % Result is in db
    disp('Computing temperature')
    T = Temperature_dev(results, conn); % Result is ITS-90 deg C
    disp('Computing conductivity')
    C = Conductivity_dev(results, conn, P.outData, T.outData); % Result mS/cm
%    disp('Computing salinity')
%     S = Salinity_dev(results, conn, P.outData, T.outData, C.outData); % Salinity, PSU, required by DO calculations
%     disp('Computing dissolved oxygen saturation concentration')
%     dOS = DOS_dev(results, conn, S.outData, T.outData); % Garcia & Gordon, ml/l
    disp('Computing dissolved oxygen voltage')
    dOV = DOV_dev(results, conn, P.outData); % Volts
    disp('Computing pH')
    PH = pH_dev(results, conn, T.outData);
    disp('Computing transmissivity')
    xmiss = Xmiss_dev(results, conn); % Beam Transmission, Chelsea/Seatech/Wetlab CStar, percent
    disp('Computing fluorescence/chlorophyl concentration')
    Chl = Fluor_dev(results, conn); % Fluorescence/Chlorophyl, Wetlab ECO-AFL/FL, mg/m^3
    disp('Computing turbidity')
    NTU = Turb_dev(results, conn); % Turbidity unit
    disp('Computing PAR')
    par = PAR_dev(results, conn); % PAR/Irradiance, Biospherical/Licor
    disp('Providing calibration data for DO')
    DO = struct('calibDate', results.Sensors.OxygenSensor.CalibrationDate,...
                'calibCoeffs', results.Sensors.OxygenSensor.CalibrationCoefficients{...
                results.Sensors.OxygenSensor.Use2007Equation+1});
  % Post-process
    for proc={'Filter', 'AlignCTD', 'CTM', 'LoopEdit', 'Derive', 'BinAverage'}
        proc = proc{:};
        reply = input(['Run ' proc '? Y/N [Y]: '], 's');
        if strcmpi(reply, 'y')
            out = [lowerCase(proc(1)) proc(2:end)];
            eval([out ' = ' proc '_dev(conn);'])
            if eval(['isprop(' out ', ''outData'')'])
                for field = eval(['fieldnames(' out '.outData)'''])
                    field = field{:};
                    eval([field 'pre' proc ' = ' field '.outData;'])
                    eval([field '.outData = ' out '.outData.' field ';'])
                end
            end
        end
    end
  % Offer to compare results to SBE Data Processing output
    reply = input('Compare to SBE results? Y/N [N]: ', 's');
    if isempty(reply), reply = 'N'; end
    
    if strcmpi(reply, 'y') % If reply was Y/y ...
        [cnvFile, cnvPath, ~] = uigetfile('*.cnv',... % Launch file chooser dialog
                                          'Choose a .cnv file');
        fid = fopen(fullfile(cnvPath, cnvFile));
        ln = fgetl(fid);
        while isempty(regexp(ln, '^(# name)'))
            ln = fgetl(fid);
        end
        
        indices = struct;
        while ~isempty(regexp(ln, '^(# name)'))
            tok = strtok(ln, '=');
            index = str2double(tok(8:end-1)) + 1;
            switch 1
                case ~isempty(regexp(ln, 'Pressure'))
                    indices.('P') = index;
                case ~isempty(regexp(ln, 'Temperature'))
                    indices.('T') = index;
                case ~isempty(regexp(ln, 'Conductivity'))
                    indices.('C') = index;
                case ~isempty(regexp(ln, 'Oxygen Voltage'))
                    indices.('dOV') = index;
                case ~isempty(regexp(ln, 'pH'))
                    indices.('PH') = index;
                case ~isempty(regexp(ln, 'Transmission'))
                    indices.('xmiss') = index;
                case ~isempty(regexp(ln, 'Fluorescence'))
                    indices.('Chl') = index;
                case ~isempty(regexpi(ln, 'Turbidity'))
                    indices.('NTU') = index;
                case ~isempty(regexp(ln, 'par'))
                    indices.('par') = index;
                case ~isempty(regexp(ln, 'Salinity'))
                    indices.('S') = index;
                case ~isempty(regexp(ln, 'Density'))
                    indices.('S') = index;
                case ~isempty(regexp(ln, 'Oxygen Saturation'))
                    indices.('dOS') = index;
                case ~isempty(regexp(ln, 'Dissolved Oxygen'))
                    indices.('DO') = index;
                case ~isempty(regexp(ln, 'Descent Rate'))
                    indices.('v') = index;
            end
            ln = fgetl(fid);
        end
        
        while isempty(ln) || ~isempty(regexp(ln, '^(*|#)'))
            pos = ftell(fid);
            ln = fgetl(fid);
        end
        [~, cnt] = sscanf(ln, '%f');
        fmt = '%f';
        for i=1:cnt-1
            fmt = [fmt ' %f'];
        end
        blocksize = 1000 * cnt;
        fseek(fid, pos, 'bof');
        SBEtmp = [];
        while ~feof(fid)
            tmp = cell2mat(textscan(fid, fmt, blocksize));
            SBEtmp = extend2Darray(SBEtmp, tmp);   
        end
        N = numel(P.outData);
        if numel(SBEtmp(:,1)) ~= N
            disp(['Sorry, FlightsDriver did not find the same number of samples in' lfcr...
                  results.dataFilename ' as were found in ' cnvFile '.  Automated comparison aborted.'])
            return
        end
        SBEdata = struct;
        for f = fields(indices)'
            f = f{:};
            SBEdata.(f) = SBEtmp(:, indices.(f));
            try
                param = eval([f '.outData']);
                SBEprec = eval([f '.SBEprec']);
            catch
                try
                    param = eval([f '.outData']);
                catch
%                     try
%                         param = eval([f '.outData']);
%                     catch
%                         try
%                             param = eval(f);
%                         catch
%                         end
%                     end
                end
            end
            switch f
                case {'P', 'PH'}
                    SBEprec = 1e-3;
%                    plotVec = SBEdata.(f) - round2(round2(param, SBEprec/10, 'fix'), SBEprec);
                    dif = SBEdata.(f) - round2(param, SBEprec);
                    if strcmp(f, 'P')
                        Title = 'Pressure';
                    else
                        Title = 'pH';
                    end
                case {'dOS', 'DO'}
                    SBEprec = 1e-5;
                    dif = SBEdata.(f) - round2(param, SBEprec);
                    if strcmp(f, 'dOS')
                        Title = 'Dissolved Oxygen Saturation Conc.';
                    else
                        Title = 'Dissolved Oxygen Concentration';
                    end
                case 'C'
                    SBEprec = 1e-6;
%                    plotVec = SBEdata.(f) - round2(round2(param, SBEprec/10, 'fix'), SBEprec);
                    dif = 10*SBEdata.(f) - round2(param, 10*SBEprec);
                    Title = 'Conductivity';
                case 'NTU'
                    sigfigs = 8;
                    SBEprec = round2(log10(abs(SBEdata.(f))), 1, 'fix') + 1 - sigfigs;
%                    SBEprec = 10^max(prec);
%                    plotVec = SBEdata.(f) - round2(round2(param, 10.^(prec-1), 'fix'), 10.^prec);
                    dif = SBEdata.(f) - round2(param, 10.^SBEprec);
                    Title = 'Turbidity';
                case 'par'
                    sigfigs = 5;
                    SBEprec = round2(log10(abs(SBEdata.(f))), 1, 'fix') + 1 - sigfigs;
%                    SBEprec = 10^max(prec);
%                    plotVec = SBEdata.(f) - round2(round2(param, 10.^(prec-1), 'fix'), 10.^prec);
                    dif = SBEdata.(f) - round2(param, 10.^SBEprec);
                    Title = 'par';
                otherwise
                    SBEprec = 1e-4;
%                     plotVec = SBEdata.(f) - round2(round2(param, SBEprec/10, ...
%                                                             'fix'), SBEprec);
                    dif = SBEdata.(f) - round2(param, SBEprec);
                    switch f
                        case 'T', Title = 'Temperature';
                        case 'xmiss', Title = 'Transmissivity';
                        case 'Chl', Title = 'Fluorescence';
                        case 'S', Title = 'Salinity';
                        case 'D', Title = 'Density';
                        case 'dOV', Title = 'Dissolved Oxygen Voltage';
                        case 'v', Title = 'Descent Rate';
                    end
            end
          % Deal w/ "bad" points
            logAbsDif = log10(abs(dif));
            logAbsDif(isinf(logAbsDif)) = ceil(min(logAbsDif(~isinf(logAbsDif)))-1);
            goodp = logAbsDif; badp = logAbsDif; reallyBadp = logAbsDif;
            if ismember(f, {'NTU', 'par'})
                good = logAbsDif < -SBEprec-1;
                reallyBad = logAbsDif > -SBEprec;
            else
                good = logAbsDif < log10(SBEprec)-1; 
                reallyBad = logAbsDif > log10(SBEprec) + 1e3*eps;
            end
            goodp(~good) = NaN; badp(good) = NaN; reallyBadp(~reallyBad) = NaN;
            figure('Name', [f ': log10(abs(SBEDP - ParseRawSBE25Data))']);
%            plot(SBEdata.(f) - param, 'k.', 'DisplayName', 'MISU unrounded') 
%            hold
            plot(goodp, '.k', 'DisplayName', ['Good Points: ' num2str(sum(good))])
            set(gca, 'Ylim', [floor(min(logAbsDif))-0.5 ceil(max(logAbsDif))+0.5],...
                     'Xlim', [1 N])
            hold;
            plot(badp, '.b', 'DisplayName', ['Bad Points: ' num2str(sum(~good))])
            plot(reallyBadp, '.r', 'DisplayName', ['Really Bad Points: ' num2str(sum(reallyBad))])            
%             plot(0*plotVec + SBEprec, 'r-', 'DisplayName', 'SBE precision')
%             plot(0*plotVec + 0.5*SBEprec, 'g-', 'DisplayName', 'Half SBE precision')
%             lwrPrec = plot(0*plotVec - SBEprec, 'r-');
%             lwrHalfPrec = plot(0*plotVec - 0.5*SBEprec, 'g-');
%             if strcmpi(f, 'DOV')
%                 ylim(SBEprec * [-1.25 1.25])
%             else
%                 ylim(SBEprec * [-1.15 1.15])
%             end
            fontSize = fix(0.8*get(gca, 'FontSize'));
            title(['Log_{10} Absolute Difference in ' Title ' Between SBEDP Software and MISU''s Conversion'...
                   lfcr 'hex: ' hexFile lfcr 'con: ' conFile lfcr 'cnv: ' ...
                   fullfile(cnvPath, cnvFile)]);%, 'Interpreter', 'none')
            set(gca, 'FontSize', fontSize)
            xlabel('Scan #')
            ylabel('log_{10}abs(SBE result - "new" result)')
            set(gcf, 'Position', [150 200 1000 600])
            lH = legend('Location', 'Best', 'Orientation', 'Horizontal');
            bestPos = get(lH, 'Position'); delete(lH); 
%            delete(lwrPrec); delete(lwrHalfPrec);
            legend('Location', bestPos, 'Orientation', 'Horizontal');
            legend('boxoff');
%             plot(0*plotVec - SBEprec, 'r-')
%             plot(0*plotVec - 0.5*SBEprec, 'g-')
        end
    else
        SBEdata = [];
    end
    fclose all;
end